NUMERICAL  INTEGRATION  EXPERIMENTS  WITH 
A  BAROTROPIC  PRIMITIVE  EQUATION  MODEL 


Jimmy  Ray  Slaughter 


VAL  POSTGRADUATE  SCHOOL 

onterey,  California 


T 


Numerical  Integration  Experiments  With 

a 
Barotropic  Primitive  Equation  Model 


by 


Jimmy  Ray  Slaughter 


Thesis  Advisor 


G.  J.  Haltiner 


March  1972 


Approved  loh.  public.  fiaXiah^',   dib&ubutlon  unlunotzd. 


Numerical  Integration  Experiments  With 

a 
Barotropic  Primitive  Equation  Model 


by 


Jimmy  Ray  Slaughter 
Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1963 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  METEOROLOGY 
from  the 

NAVAL  POSTGRADUATE  SCHOOL 
March  1972 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CALIFORNIA   93943-5002 


ABSTRACT 

Experiments  are  carried  out  with  various  time  and  space 
differencing  schemes  applied  to  the  barotropic  primitive 
equations  using  both  real  data  and  a  particular  stream  func- 
tion which  is  an  analytic  solution  to  the  nondivergent  baro- 
tropic vorticity  equation.   With  both  types  of  data,  there 
were  some  significant  differences  in  the  forecasts  produced 
by  the  various  schemes.   Replacement  of  the  widely  used  leap 
frog  (centered)  scheme  by  others  which  eliminate  or  lessen 
some  of  its  inherent  errors  at  the  expense  of  more  computer 
time  or  storage  appears  to  be  justified  at  such  time  when 
computer  capacity  no  longer  restricts  operational  use  of 
these  more  time  consuming  schemes. 
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I.   INTRODUCTION 

In  selecting  a  particular  differencing  scheme  for  appli- 
cation to  numerical  weather  prediction,  a  major  considera- 
tion is  computer  time  and  storage  requirements.   Therefore, 
the  selected  scheme  is  a  compromise  between  verification  ac- 
curacy and  computer  limitations.   At  present,  the  leap  frog 
(centered)  scheme  is  widely  used  in  numerical  weather  predic- 
tion and  will  be  used  as  a  standard  for  purposes  of  compari- 
son.  Present  computer  limitations  and  methods  of  verifica- 
tion tend  to  preclude  seeking  more  accurate  solutions  based 
on  mathematical  considerations  alone.   The  purpose  of  this 
study  is  to  look  beyond  present  computer  capabilities  and 
examine  several  other  differencing  schemes  which  reduce  or 
eliminate  some  of  the  inaccuracies  in  the  leap  frog  scheme. 
This  scheme  has  three  time  levels  and,  therefore,  a  spurious 
computational  mode  arises  from  the  finite  differencing 
scheme.   This  problem  can  be  eliminated  by  the  application 
of  two  level  time  schemes  such  as  the  Euler  backward  or 
Runge  Kutta  methods  which  have  on  computational  mode.   A 
fourth  order  space  differencing  scheme  is  applied  to  reduce 
truncation  error  and  improve  phase  speeds.   Should  any  of 
the  schemes  produce  a  significant  change  in  forecast  rela- 
tive to  the  leap  frog  method,  it  must  be  determined  if  the 
difference  represents  a  sufficient  increase  in  accuracy  to 
warrant  the  added  computer  time  and  storage  requirements. 


II .   PROGNOSTIC  EQUATIONS 

The  equations  of  motion  and  the  continuity  equation  for 
a  barotropic  fluid  in  the  momentum  flux  form  are 


9(uh)  2  9  (uuh)^  2    9  (uvh)  _  ,  v9h 

fr-r +  m   -r—  —  +  m   "^ —  fvh  =  -  mgh-x — 

dt  dx  m        ay  m      dx 

9(vh)  .  2  8  (vuh)J  2  8  (vvh)±  ,  ,  ,9h 

~ +  m   -r—  —  +  m   tt—  +  fuh  =  -  mgh-r— 

9  t  9  x            9  y  m  —  9  y 

9h     2  9  (uh)    9  (vh),  n 

•z—   +  m  [it—   +  -z —  J  -  0 

9t  9xm  dy  m 


(1) 
(2) 
(3) 


where  h  is  the  depth  of  the  fluid  (or  height  of  a  pressure 
surface  in  the  atmosphere),  g  is  the  force  of  gravity,  f  is 
the  coriolis  parameter  and  u  and  v  are  horizontal  velocity 
components . 

A.   SPACE  DIFFERENCING 

1 .   Second  Order  Scheme 

With  u,  v,  and  h  given  at  each  grid  point,  the 
differential  equations  may  be  approximated  by 


uh   -   .   uh  .  . 

0  (uh)  .  .        2//jr/  _l       w    1  +  1>3.      1J\ 

it  =  -m  /4df (u. .,   .+  u..)( 'J  +  J) 

m         m 


9t 


3  -  -m  /Ad[(u.+1>.+  u..)(; 


(4) 


uh 


uh .  - 

-  (u.  .   .+  u. .)( X  1,3+  1J 

l-l, j    13   m         m 


) 


vh .   .  .  n   vh  .  . 

+  (u.  .._+  *       )t-J-*+1+   — 1J) 

1,3+1    13   m         m 

vh .  .  '    vh  .  . 

-  (u.    n  +  u..)( 1,J  i+  1J)]  +  f..(vh).. 

1,3-1    13   m         m  13    '13 

-  mgh.  .  (h.  ,  .   -  h.  .   .)/2d  =  F.  . 

13   1+1,3    i-l»3  13 

where  3  increases  with  x  and  i  decreases  with  y.   (See  Fig. 
1)   Similar  expressions  are  used  for  the  other  two  equa- 
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tions.   This  method  of  space  differencing,  devised  by  A.  Ara^ 
kawa  (see,  for  example,  [1]),  prevents  nonlinear  instability 
in  the  continuous  time  case  (sometimes  referred  to  as  the 
semi-discrete  equations). 

2 .   Fourth  Order  Scheme 

Fourth  order  space  differencing  which  takes  into 
account  a  greater  number  of  grid  points  and  reduces  trunca- 
tion error  approximates  uh  tendencies  by 

2 


f^l.3   -  */3(,1J) 

uh 


-H- 


[<u  +    u..)(— i+2'j        (5) 

8d  i+2, j  xj       m 


m 


li  JJ 


m 


)-(u,     ,     4+    u.4)< 


+  (u 


i-2,j  ij       m 

vh .  vh . . 


uh .     n     .       uh .  . 
1-2, J+   : 


m 


) 


i,j+2         ij       m 


m 


vh  .  _       vh .  . 

-<u.  +    u..)(— X»J     Z+   — 1J)]    +    f       (vh).. 

i»J-2  ij       m  m  ij  xj 

-mgh.  .  (h,  ,_     .     -    h.  ,)/4dl 

xj       x  +  2,j  i-2, j  J 

where  F.  .  is  the  second  order  space  difference. 
i,J 

B.   TIME  DIFFERENCING 

U  represents  uh ,  vh ,  or  h;  t  is  time;  n  is  time  step 

r)  TT 

number  and  F  is  -^ —  in  the  following  time  differencing 

o  t 

s chemes  : 

1 .  Leap  Frog 

Un+1  _  u11-1  =  2AtFn 

2 .  Runge  Kutta 

AUi  =  At  F[Un,  nAt] 


AU 


2  =  At  F[Un+  AU1/2,  (n+j)At] 
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AU 
AU3  =  At  F[Un  +  ——    ,  (n  +  p    At] 


AU4  =  At  F[U   +  AU   ,  (n  +  1)  At] 


Un+1  =  Un  +  rt^U-  +  2AU0  +  2AU_  +  AU.] 

O      1  l  J         4 


3.   Euler  Backward 


U*  -  Un   =   At  Fn 
Un+1-Un   =   At  F* 


Denotes  trial  values. 
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Ill .   BOUNDARY  CONDITIONS 

A.   A  PARTICULAR  STREAM  FUNCTION 

The  boundary  conditions  in  this  case  are  applied  to  an 
initial  stream  function  which  is  an  analytic  solution  to  the 
nondivergent  barotropic  vorticity  equation  and  after  each 
time  step  during  the  forecast. 

1 .  North-South   Boundary  Conditions 

For  the  63x63  grid  a  wall  exists  between  y=l  and 
y=2  and  also  between  y=63  and  y=62  (See  Fig.  1).   No  flux 
is  permitted  across  this  rigid  boundary.   This  is  accom- 
plished by  setting  v   ,  =  -v   ,_  and  v    =  -v   0  thereby 

x,63     x,62       x,l     x,2 

making  the  average  v    =0  alone  the  wall.   The  other  para- 

x  ,  y 

meters  are  set  equal  across  the  wall.   Summarizing,  the 

boundary  conditions  are 

North  South 

v   0  =  -  v   ,,  V   n  =  -  V   «, 

x,  63       x  ,  62  x,l       x,2 

U     (LI 

x ,  63  =    u   , 0  un=u0 

x,  62  x  ,  1       x , 2 

*x,63  =    *x,62  *x,l  =    *x,2 

where  (j)   is  the  geopo tent ial . 

2 .  East-West  Boundary  Conditions 

The  East-West  boundary  conditions  are  cyclic  as 
suggested  by  the  function  plotted  in  Fig.  1.   Summarizing, 
the  boundary  conditions  are 


13 


East  West 

l,y     62, y  63, y     2,y 

where  A  represents  parameters  v     ,  u     and  d> 

x,y     x,y      Yx,y 

3 .   Variation  for  Fourth  Order  Space  Difference 

In  order  to  apply  fourth  order  space  differencing 
and  still  maintain  the  boundary  conditions  of  the  original 
grid,  it  was  convenient  to  create  an  additional  outside  row 
of  points  and  relabel  as  a  65x65  grid.  (See  Fig.  2)   The 
boundary  conditions  are  summarized  as  follows: 

North  South 

x,65     x,62         x,l      x,4 

A     represents  u    , 

x>y     j  a       x,y 
v   ,.  =-v   ,_        v   „  =  -v   0         J       and  m 
x,64     x,63         x,2      x,3  x,y 

x,65    x,62         x,l      x,4 

x ,  54    x,63         x,2      x,3 

East  West 

A      =  A  A     =   A 

l,y      62, y         65, y     4,y      A     represents  u    , 

x  >  y  x  >  y 

A0     =  A-„  A,.      A_        v     ,  and  (f> 

2,y      63, y         64, y     3,y       x,y        Yx,y 

B,   ACTUAL  500  MILLIBAR  DATA 

When  the  equations  were  integrated  using  real  data,  the 

"spongy"  boundary  conditions  of  the  FNWC  model  were  applied 

After  each  time  step,  the  prognostic  results  were  modified 

by  setting  t  f .        , 

uh  =  (1-k  )uh_  +  k   uh  b 

vh  =  (l-k')vh  +  k'  vh0 

h  =  (l-k')h   +  k'  h0 
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where  the   0  subscript  indicates  initial  values.   Below  lat- 
itude four  degrees  north  set  k'=4.   Above  17  degrees  north 
set  k'=0.   Between  four  degrees  north  and  17  degrees  north, 
there  is  linear  variation  in  k '  from  zero  to  one.   This  re- 
stores the  variables  to  their  initial  values  south  of  four 
degrees  north,  while  north  of  17  degrees  north,  full  varia- 
bility is  permitted.   There  is  no  restoration  north  of  17 
degrees  north.   Between,  there  is  partial  restoration. 
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^_     y  =  £3    •> 

/  / 1 1  /it/ 1 it/ zrnxz&L   *wfl[*  "^xxaxxxzznxzn  zzzzzzzzzzzgafe 


y»tf2  ^ 


r 


X«l 


*2 


y-2. 


iu/i/i//i//ii ///tuff**  ^JJft    tg./i  i  mui/  uiiiinii  wtr/u,// 


Illustrating  the  63x63  square  grid  and 
boundaries  as  applied  to  a  particular 
stream  function.   The  dimens ionless  sym- 
bols x  and  y  are  used  for  notation  con- 
venience . 

Figure  1 
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~-   H=GS 


/VORTH 


rn/i  >//>  jf/t  j  j/j  ///jrfZ>    WfiU     ?&/)}  //  /  /  j/f  sit  >  >  //  /  fti-r? 


»j--63 


y--62 


/ /JtJj/7 


X=l        X«2        *3       Xsif 

1/1  i   f 


i   \/\ 


X-S2      X'G3     Xs  C¥     X'fer 


!/ 


I 


y=* 


y-3 


$ 


gzzzzzgzzzzzzzz  zzzzzzzzzz^c  5$JJ  ^zazzzzzzzzzzzzzzzzzzzzzzzzgzzzzzzz 

yr2. 


J 


Illustrating  the  65x65  expansion  of  the  grid  in 
Fig.  1  in  order  to  apply  fourth  order  space  dif- 
ferencing over  the  same  area.   The  dimens ionles s 
symbols  x  and  y  are  used  for  notation  convenience 

Figure  2 
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IV.   GRID 

A.  A  PARTICULAR  STREAM  FUNCTION 

A  63x63  square  grid  as  in  Fig.  1  and  as  modified  in 
Fig.  2  for  fourth  order  space  differencing   is  used  to 
minimize  the  number  of  changes  required  to  adapt  programs 
to  the  FNWC  63x63  hemispheric  grid.   The  boundaries  in 
latitude  extend  from  30  to  90  degrees.   Grid  length  is 
determined  by  dividing  the  geographic  distance  represented 
by  this  latitude  difference  into  62  equal  parts. 

B.  FNWC  63x63  HEMISPHERIC 

In  conjunction  with  the  actual  500  millibar  data,  the 
FNWC  63x63  hemispheric  grid  is  used.   This  is  a  square  re- 
presentation of  a  hemispheric  polar  stereographic  grid  in 
which  the  grid  distance  d  is  determined  by 

,     d1    , .  .  ,1  +  sin  60  * 

d  =  "m-  =  d  ' (1  +  sin  Y   } 

where  m  is  the  map  factor,  d'  is  the  grid  distance  at  60 

latitude  (381  kilometers),  and  sin  y  is  defined  by 

973.71  -  R2    ,     D2    ,_  _  .2,  ,_  _  ,2 

sxn  y  =   2   wnere  R   =  (I-I  )  +  (J-J  ) 

973.71  +  R  P         P 

I, J  =  arbitrary  grid  point. 

I  ,J  =  grid  point  for  the  north  pole  =  32,32. 
P   P 
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V.   INITIAL  CONDITIONS 

A.   A  PARTICULAR  STREAM  FUNCTION 

In  order  to  obtain  a  smooth  initial  stream  function 
field,  an  analytic  solution  to  the  nondivergent  barotropic 
vorticity  equation  is  used,  namely 

]b    -    lb       sin  - —  (y  -  -r)  cos  - —  (x-ct) 
°       L„        I  L 

y  x 

where  c  is  the  phase  speed  and  the  wavelength  in  the  x  dir- 
ection and  the  half  wavelength  in  the  y  direction  are  equal 
to  61d  (L„  =  L„  =  61d) .   The  initial  stream  function  is 

a  y 

i  i  T       /  d>  2tt 

y    =    i>      sin   - —    (y    -    -r-)     cos    - —   x    . 

O  Li  Z  L„ 

y  X 

Here,  \p       is  determined  by  using  the  relation 

v  L 
I    z  o  X 

^o      2tt 

where  v   is  the  initial  specified  amplitude  of  the  meri- 
o 

dional  wind  speed  derived  from 

drp  2tt  ,    r  .   IT   ,    dN       2  tt   , 

V  =   9^   =  L~  ^o  [sin  L~  (y"  2>  COS  Z~   X]* 
x  y  x 

Using  the  initial  field  thus  obtained,  the  linear  balance 

equation 

V2<f>    =    fV2ijj   +   Vf    •    ViJ; 

which  in  finite  difference  form  may  be  written 

?2<f> .  .=    f  .  .    V2iK.    +t[^    f    v    ^+    v    f    V    i>]    =    F.  . 

1J  1J  Y1J  4         X  X  y  yyj  ij 
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where 

*2<    >  -  <  >i+i(j+  <  >t-i.i+  <  >i,j+i+  <  >i.j-i-  4<  >u 

\-   <  >i+i.3  -  <  h-l,i   ««=• 

The  preceding  equation  may  now  be  solved  as  a  Poisson 

equation  for  the  geopotential  field.   Thus, 

V2<K  .  =  F.  . 
ij     ij 

The  initial  guess  for  the  relaxation  procedure  is 

6   =  f  \b  where  f  is  a  mean  coriolis  parameter. 

To      r  o  r 

Subsequent  "guesses"  are  made  according  to  the  relation 


ij 

where 

R.(n)    =    a(?2<J>.fn)  -  F.  .) 
ij  ij       ij 

and  a  is  a  relaxation  coefficient. 

Using  h  =  4>/g,  computations  for  the  initial  height  field 
are  completed.   Additionally,  a  constant  height  is  added 
to  each  grid  point  value  in  order  to  eliminate  negative 
heights,  which  would  otherwise  cause  computational  insta- 
bility. (See  Fig.  3) 

B.   ACTUAL  500  MILLIBAR  DATA 

The  initial  height  field  is  obtained  from  the  FNWC  500 
millibar  analysis.   The  linear  balance  equation  is  solved 
for  the  stream  function  as  a  Poisson  equation  using  ip=  <J)/f 
as  a  first  guess . 
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Differencing 
Scheme 

Leap  Frog 

Leap  Frog  Fourth 
Order  Space 

Runge  Kutta 

Euler  Backward 


RMSE  For  Position 
(Nautical  Miles) 


RMSE  For  Central 
He i  ght  (Me t  e  r  s ) 


97 

30 

76 

31 

93 

31 

88 

25 

Table  of  root  mean  square  errors  in  position  and  central 
height  of  the  three  low  centers  at  24,  48,  and  72  hours 
as  forecast  by  various  differencing  schemes  (on  actual  data) 
in  comparison  to  the  corresponding  analysis. 

TABLE  1 


Latitude  =  90 


{  Goo   580     WO  tOO 


520  600    J-j 


Latitude  =  30 


Initial  height  field  as  produced  by 
cular  analytic  stream  function  on  a 
grid.  Only  the  central  contours  of 
features  are  illustrated. 

Figure  3 


a  par t i- 

63x63 
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VI.   RESULTS 

The  results  shown  in  the  following  figures  were  obtained 
from  the  computer  subroutine  mapping  of  the  63x63  height 
field  obtained  by  the  various  differencing  schemes  using  10 
minute  time  steps.   One  grid  squares  are  used  in  the  figures 
simply  as  a  convenience  to  provide  a  ready  reference  scale 
to  compare  positions  of  low  centers. 

A.  PARTICULAR  STREAM  FUNCTION 

Results  for  the  forecasts  made  using  a  particular 
stream  function  for  the  initial  values  are  summarized  in 
Fig.  4  for  central  height  and  relative  position.   The  fore- 
cast results  for  the  individual  differencing  schemes  are 
presented  in  greater  detail  in  Fig.  8  through  Fig.  12. 

B.  ACTUAL  500  MILLIBAR  DATA 

The  results  of  forecasts  made  on  three  low  systems  from 
the  500  millibar  FNWC  analysis  are  summarized  in  Fig.  5, 
Fig.  6  and  Fig.  7.   Although  only  low  systems  are  illustra- 
ted in  this  study,  similar  results  were  obtained  from  the 
movement  and  intensification  of  highs.   The  time  frame  of 
this  study  begins  with  the  1200Z  March  24,  1966  500  milli- 
bar heights  and  terminates  with  the  1200Z  March  27,  1966 
forecast.   Low  number  one  (See  Fig.  5)  was  located  over 
central  Siberia.   Low  number  two  (See  Fig.  6)  and  low  num- 
ber three  (See  Fig,  7)  were  located  over  Southwestern 
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Canada  and  Labrador  respectively.   Root  mean  square  error 

computations  for  position  and  central  height  as  forecast 

by  each  differencing  scheme  when  compared  to  the  correspond- 
ing analysis  are  contained  in  Table  1. 

C.   COMPUTER  REQUIREMENTS 

Although  peak  program  efficiency  was  not  attempted  dur- 
ing the  application  of  the  various  numerical  schemes,  con- 
sistency in  application  was  carefully  maintained.   Using 
the  leap  frog  scheme  as  a  standard,  the  leap  frog  fourth 
order  space  differencing  and  Euler  backward  schemes  require 
approximately  a  35  percent  increase  in  run  time.   The  Runge 
Kutta  scheme  requires  approximately  a  170  percent  increase 
in  run  time  and  a  25  percent  increase  in  computer  storage. 
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Figure  7 
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VII.   CONCLUSIONS 

The  various  differencing  schemes  were  applied  to  the 
smooth  analytic  solution  in  order  that  differences  in  fore- 
casts produced  by  these  schemes  would  be  more  readily  dis- 
tinguishable than  with  actual  data.   However,  the  choice 
as  to  resolution  was  probably  too  high  in  this  study  (61 
grid  lengths  per  wavelength)  and  shorter  wavelengths  are 
yet  to  be  tested.   Differences  produced  in  location  and 
central  values  were  not  conclusive  when  compared  to  the 
results  produced  using  actual  data.   A  reduction  in  resolu- 
tion for  the  analytic  initial  field  (approximately  six  grid 
lengths  per  wavelength)  would  possibly  produce  a  more  mean- 
ingful separation  in  the  accuracy  produced  by  the  various 
schemes . 

The  negligible  increase  in  accuracy  (four  percent  based 
on  a  root  mean  square  error  comparison)  of  the  Runge  Kutta 
scheme  over  the  leap  frog  scheme  does  not  justify  the  ex- 
tremely large  increase  in  computer  time  and  moderate  in- 
crease in  computer  storage.   Both  the  leap  frog  fourth 
order  space  and  Euler  backward  differencing  schemes  show 
promise  for  improved  accuracy  with  a  moderate  increase  in 
computer  run  time  (approximately  35  percent) .   The  leap 
frog  fourth  order  space  differencing  scheme  gave  a  24  per- 
cent increase  in  position  accuracy.   The  Euler  backward 
•scheme  produced  a  nine  percent  increase  in  position  accur- 
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acy  and  a  17  percent  increase  in  central  height.   With  both 
the  Euler  backward  and  the  leap  frog  fourth  order  space 
schemes,  the  35  percent  increase  in  computer  time  could  be 
reduced  somewhat  by  more  efficient  programs.   The  small 
statistical  sample  of  pressure  systems  examined  in  this 
study  does  not  permit  firm  conclusions  to  be  drawn  regarding 
the  various  schemes  tested.   However,  there  is  reason  to 
believe  that  the  fourth  order  space  differencing  of  the 
nonlinear  advection  terms  will  reduce  the  phase  error  in  the 
prediction  of  pressure  systems  and  this  method  should  be 
tested  further.   Williams  [2]  has  shown  that  the  pressure 
force  need  only  be  differenced  with  second  order  accuracy 
which  permits  a  somewhat  larger  time  step  without  computa- 
tional instability  than  when  fourth  order  differences  are 
used  throughout. 
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